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We report a study of viscosity by the method of Helfand moment in systems with periodic bound- 
. ary conditions. We propose a new definition of Helfand moment which takes into account the 

minimum image convention used in molecular dynamics with periodic boundary conditions. Our 
Helfand-moment method is equivalent to the method based on the Green-Kubo formula and is not 
affected by ambiguities due to the periodic boundary conditions. Moreover, in hard-ball systems, 
our method is equivalent to the one developed by B. J. Alder, D. M. Gass, and T. E. Wainwright 
[J. Chem. Phys. 53, 3813 (1970)]. We apply and verify our method in a fluid composed of N > 2 
hard disks in elastic collisions. We show that the viscosity coefficients already take values in good 
agreement with Enskog's theory for N — 2 hard disks in a hexagonal geometry. 
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d . I. INTRODUCTION 

(— > 

Viscosity is the fundamental mechanism of dissipation of momentum in a fluid. Viscosity is described at the 
macroscopic level by the Navier-Stokes equations which are the equations of balance of momentum in a fluid. At 
the microscopic level, viscosity arises because of a transfer of momentum between fluid layers moving at different 
' velocities as already explained by Maxwell thanks to kinetic theory. 

In the fifties, Green, Kubo, Mori, and others provided an explanation of all the transport properties in terms of 
time-dependent statistical correlations of microscopic currents associated with each transport property 0, 0, El 0] • 
O , They showed that the transport coefficients are given as the time integrals of the time autocorrelation functions of 
the microscopic currents, yielding the famous Green-Kubo formulas. Thereafter, Helfand showed in the early sixties 
y—i \ that the transport coefficients can be expressed by Einstein-like formulas in terms of moments - the so-called Helfand 
^ ■ moments - which are the time integrals of the microscopic currents p| . 

These new methods by Green, Kubo, Mori, Helfand, and others have been applied to the computation of transport 
CM properties by molecular-dynamics simulations, in particular, by Alder et al. [6j. In molecular-dynamics simulations, 
the system is necessarily composed of a finite number of particles which are usually moving in a domain defined with 
periodic boundary conditions in order to simulate the bulk properties. The periodic boundary conditions (p.b.c.) 
■ usually considered in molecular dynamics are based on the so-called minimum image convention according to which 
interaction should occur between pairs of particles separated by the minimum distance among the infinitely many 
"jS i images of the particles allowed by the p.b.c. In molecular-dynamics simulations, the minimum image convention 
plays a fundamental role to define the microscopic current entering the Green-Kubo formula. 

We may wonder if the Helfand-moment method could be applied to molecular dynamics simulations with p.b.c. 
The advantage of the Helfand-moment method is that it expresses the transport coefficients by Einstein-like formulas, 
directly showing their positivity. Moreover, this method is very efficient because it is based on a straightforward 
accumulation which is numerically robust. Actually, it is a Helfand-moment method which has been numerically 
q | implemented by Alder et al. for viscosity in hard-ball fluids Several other implementations of the Helfand- 
moment method have been considered and discussed in the literature However, the implementation of this 
. £h ' method for systems subject to p.b.c. other than hard-ball fluids seems to remain ambiguous as reported by Erpenbeck 
^ ! in Ref. @. 

The purpose of the present paper is to propose a Helfand-moment method which is appropriate for molecular 
dynamics simulations with p.b.c. and which is strictly equivalent to calculations with the Green-Kubo formula. For 
this purpose, we show the need to take into account the minimum image convention. In this way, we are able to 
obtain a Helfand moment giving viscosity thanks to an Einstein- like formula in molecular dynamics with p.b.c. The 
so-obtained value of viscosity is in full agreement with the value of the Green-Kubo formula and also with the value 
obtained by Alder et al. @- 

Our method is applied to the hard-disk fluid. We study in detail the simple model composed of two 
hard disks in elastic collisions in a domain defined by p.b.c. Due to the defocusing character of the disks, 
this model is chaotic. Bunimovich and Spohn have demonstrated that the viscosity already exists in this 
system with only two particles 0. The model they studied is defined with p.b.c. in a square geom- 
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etry. It presents a fluid and a solid phases which are separated by a phase transition. The problems 
presented by the model in a square geometry are that: (i) the viscosity exists only in the solid phase; 
(ii) the viscosity tensor which is of fourth order is anisotropic on a square lattice. In the present work, we solve 
these problems by considering a hexagonal geometry. Indeed, in the hexagonal geometry, the fourth-order viscosity 
tensor is isotropic and we can proof the existence of viscosity already in the fluid phase. 

Furthermore, we apply our method to systems containing more and more hard disks. We show that the values of 
the shear viscosity obtained by our Helfand-moment method are in good agreement with Enskog's theory, already for 
the fluid of two-hard disks. 

The paper is organized as follows. In Sec. [DJ we derive our expression of the Helfand moment for the viscosity 
tensor to be of application in molecular-dynamics simulations with p.b.c. In Sec. 11111 we describe the model of two 
hard disks in the hexagonal and square geometries. In Sec. IIVI we study different properties of the model like the 
mean free path and the hydrostatic pressure, in particular, across the fluid-solid phase transition. In Sec. the 
Helfand-moment method is applied to the two-disk model to calculate the shear and bulk viscosities. We show how 
the fluid-solid phase transition affects the viscosities in this model. In Sec. IVII we extend our results to systems 
with N — 4,8, 12, ...,40 hard disks. We show that the shear viscosity already takes a value in good agreement with 
Enskog's theory in the two-hard-disk system. Our results are discussed and conclusions are drawn in Sec. IVIII 



II. HYDRODYNAMICS, HELFAND MOMENT, AND VISCOSITY 

A. Viscosity and hydrodynamics 

The hydrodynamic theory provides us with the equations of motion for the conserved quantities in a fluid. In 
particular, the local conservation of momentum is expressed by the well-known Navier-Stokes equations |ll|: 

dpvj _ dTljj . , 

dt ~ dTj ' 1 J 

where 



II., = pViVj + P 8ij - o\, , (2) 

is the momentum flux density tensor, pvi is the momentum density, P the hydrostatic pressure, and a'^ the viscosity 
stress tensor. This last tensor takes into account the internal friction occurring in a fluid when different parts of the 
fluid move with different velocities. Therefore, a'^ has to be proportional to the space derivatives of the velocities: 

o-ij = T)ij,kl -q^ , (3) 

where rjij^ki is the viscosity tensor. 

For isotropic systems, the theory of Cartesian tensors shows that the basic isotropic tensor is the Kronecker tensor 
Sij and that all the isotropic tensors of even order can be written like a sum of products of tensors 5ij |12| : 



Vij,ki = a S i3 Sm + b 8ik Sji + c Sjk Su , (4) 

where a, b and c are scalars. Since the viscosity stress tensor is symmetric a'^ = cr^, only two of these coefficients are 
independent because b = c. After a rearrangement, we have: 

, / dvi dvj 2 c dv[\ A c dvi 



^=H^ + ^-^J +C ^' (5) 

for a d-dimensional system. The coefficients r/ = b and £ = a + (2/d)b are respectively the shear and bulk viscosities 
and they can be expressed in terms of the elements of the fourth-order viscosity tensor as: 

I C ~j T]xx,xx ~} T]xx,yy • ^ ^ 
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B. The Green-Kubo formula in molecular dynamics with p.b.c. 

Several techniques have been developed during the last century to evaluate the transport coefficients. One of the 
most important methods was established by Green, Kubo, and Mori 0, 0, 01- It consists in having a relation 
between each transport coefficient and the autocorrelation function of the associated flux or microscopic current. In 
our case (see Appendix lA"|l . we have: 



W = | / [(J«(0) Jki{t)) - {JiMJki)] dt , (7) 



with the microscopic current: 



N N 

Jij = E Vj ^r + 2 £ Fi(r - - n) ~ r ^ ' (8) 

a—l 

p a and r a being the momentum and the position of the a th particle, while F(r a — r ) is the force between particles 
a and b. In Eq. J7J the average (•) is performed with respect to the equilibrium state. We notice that, for the 
microcanonical state, 

(see Appendix [XJ. 

A very important point is that, in a system with p.b.c. as considered in molecular-dynamics simulations, the 
difference of positions r a — r;, must satisfy the minimum image convention that 

\r a j-r bj \<^, for j = l,...,d, (10) 

for a cubic geometry. More generally, the difference of positions must remain within a unit cell of the Bravais lattice 
used to define the p.b.c. With p.b.c, there is indeed an infinite lattice of images of each particle. All these images 
move in parallel. If the force has a finite range the particle a interacts only with the particles b within its interaction 
range The force field F(r) has a finite range of interaction beyond which it vanishes. The interaction range is 
supposed to be smaller than the size L of the box containing all the particles. It is important to notice that we do not 
suppose here that the force field is periodic. In order to define a dynamics which is periodic in the box of size L the 
positions should jump in order to satisfy the minimum image convention. As a consequence of this assumption, the 
positions and momenta used to calculate the viscosity by the Green-Kubo method actually obey modified Newton's 
equations 

£_£ +5>f > ,(«-«.), 

6(/a) 

where Ari^ is the jump of the particle a at time t s in order to satisfy the minimum image convention. Moreover, 
we impose that the particle No. 1 does not jump. To satisfy these conditions, the jumps at the time t s when 
\r a j(t s ) - r bj (t s )\ = L/2 can be given by 



for a < b 



Ar# = , 
Arlf=eL, 

Ar ( c f = , for c ^ a, b , 
Ar^ = , for k ^ j and Vd 



(12) 



with e = sgn[p a j(t s ) —pbj(t s )]. The modified Newton equations (|11|) define a dynamics which is periodic on the torus 
of the relative coordinates r a — ri because the jumps of the relative coordinates are vectors of the Bravais cubic 
lattice: Ar a ^ — Ar^ = 0, ±L, while the momenta p a remain functions of the time without singularities worst than 
discontinuities. We notice that modified Newton's equations 1)110 conserve energy, total momentum and preserve 
phase-space volumes (Liouville's theorem). 
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C. Helfand moment for molecular dynamics with p.b.c. 

In the sixties, Helfand has derived quantities associated with the different transport processes, in particular for the 
viscosities These new quantities Gij(t) are such that we can obtain an Einstein- like relation for each transport 
coefficient. For the shear viscosity coefficient, we have: 

V = t Hm ([G xy (t) - G xy (0)f) ■ (13) 
More generally, we can define such a relation for each element of the viscosity tensor: 

^ = & 4v [{Gi] {t)Gki {t)) - {Gt] {t)) {Gki m (i4) 

if we take (7^(0) = 0. The Helfand moment Gy (i) is defined as the integral of the microscopic current appearing in 
the Green-Kubo relation: 

Gy(t)=G«(0) + I M r ) dT ■ ( 15 ) 
Jo 

As a consequence of the definition (|f 5JI . the Einstein-Helfand formula Ijf 4|l is equivalent to the Green-Kubo formula 
, as proved in Appendix [BJ In a system of N particles on a torus and satisfying the minimum image convention, 
we can integrate the current © with modified Newton's equations l|II|) to get: 

N N 

Gij(t) = $>a«(i) r aj (t) - Ar# 9(t - t s ) , (16) 

a—1 o— 1 s 

where £^(0) = 0, p^ — p a i(ts) — t s ) is the Heaviside step function at the time t s of the jump s: 

The expression l|I6|) which we propose in the present paper can be used to obtain the viscosity coefficients thanks to 
the Einstein-like formulas (114ft in a molecular dynamics defined on the torus. We emphasize that the expression l)16|l 
may apply to systems of particles interacting with a smooth potential under the condition that the range is finite, 
or to systems of hard balls in elastic collisions. We show in Appendix [C] that the hydrostatic pressure can also be 
written in terms of the Helfand moment (|16J) . 

Our Helfand- moment method has several theoretical and numerical advantages: (i) It is strictly equivalent to 
the Green-Kubo method, (ii) The Einstein-like formula l|13|l or (f 1 4f) directly show the positivity of the viscosity 
coefficient or viscosity tensor because t, /3, and V are positive. Moreover, the Helfand moments directly obey central 
limit theorems, expressing the Gaussian character of the dynamical fluctuations in systems with finite viscosity, 
(iii) Thanks to our expression (|16|1 of the Helfand moment, the viscosity coefficients are given by a straightforward 
accumulation over the successive jumps s. For a given system with N particles, numerical convergence can be reached 
in the limit of an arbitrarily large number of jumps s, under conditions of existence of the viscosity coefficients. 

By defining the Helfand moment as the integral (|15|l of the microscopic current for a system with minimum image 
convention, we obtain the expression (|16|l which can be used to directly calculate AGy (t) = Gij (t) — Gij (0) for the 
Einstein-Helfand relation, remaining consistent with the requirements imposed by the periodic boundary conditions 
and with the Green-Kubo formula for a system satisfying the minimum image convention. 

D. Comparison with other methods 

In the seventies, Alder et al. 0] calculated the viscosity coefficients of hard-ball systems with Einstein-like formulas 
based on expressions for Helfand moments which are specific to hard-ball systems. Instead of adding a new quantity 
to the Helfand moment at each passage through the boundaries of the minimum image convention as in Eq. Ultjfl > 
their expression takes into account only the elastic collisions between the hard balls. The Helfand moment can be 
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obtained by direct integration of the microscopic current according to Eq. (|15(l with Gy(0) = 0: 

Gij(t) = fdrJijir) (18) 
Jo 



t 

dr 



o 



N 1 

E + 2 E F *( r " " rfc) ~ rb ^ 

o=l a^b 



(19) 
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Between the collisions the trajectory is a straight line and the particle velocities change only at each collision. There- 
fore, the first term in the integral, the kinetic term, is constant during two successive collisions and changes only at 
the collisions. The second term, the potential term, vanishes between two successive collisions and contributes only 
at collisions. Indeed, for a hard-ball potential, the forces between the particles a and b colliding at the time t c of the 
collision c can be written in terms of the change Ap„ = p a (t c + e) — Pa{t c — e) of momentum of the particle a at the 
collision c as 

-Ap^ S(t t c ) , 
-Api c) 6(t - t c ) , 

(c) (c) 

for t c — e < t < t c + e, because Ap^ = — Apa . The forces with the other particles which are not engaged in the 
collision vanish. Therefore, we obtain: 

GiAt)= E (E^| /^v + E^ff (21) 

(c-l,c) 

where, in the first term, At c _i. c is the time of flight between the collisions c— 1 and c during which the momenta remain 
constant and, in the second term, a and b denote the particles interacting at the collision c and r^- = r a j{t c ) — ^bj{tc) ■ 
The first sum runs over the intercollisional free flights (c — 1, c) between the initial time t = and the current time t, 
while the second sum runs over the collisions occurring between the time t = and t. If C denotes the last collision 
before the current time t, we notice that the last term of the first sum is Atc,c+i =t — tc- Hence, if we differentiate 
Eq. (|21|l with respect to time and use l|2(JII we recover the microscopic current (JSJ. Therefore, the expression (|21ll is 
equivalent to our expression l|16l) in the case of hard-ball systems. However, our expression (|16|l extends to systems 
with a smooth interaction potential. 

A comment is here in order about another method which has been considered and discussed in the literature 
0, El ■ This other method implements an expression printed in the middle of a presentation given in Ref . |13| for 
the calculation of shear viscosity with the Helfand-moment method. This expression differs from the Helfand moment 
by the mere exchange of a square and a sum over the particles. The equivalence of the expression in Ref. ^| with 
the Helfand moment depends on the vanishing of some cross terms as pointed out in Ref. y\ . Numerical evidence has 
been obtained in Refs. |8|, |9( that the expression in Ref. |13( is in general not valid to calculate the shear viscosity. 
We notice that both the original Helfand moment and the expression in Ref. do not strictly apply to systems 
subject to p. b. c. (see discussions in Refs. @, @)- This problem is solved by the expression (|2T)l of Ref. [(| in the 
case of hard-ball fluids and by our expression (|Pd|) in the general case. 



E. Symmetry considerations in two-dimensional systems 



By symmetry, most of the elements of the viscosity tensor are either equal or vanish. First, we have: 

Vij,kl — Vkl,ij — Vji,kl = Vij,lk , (22) 

because of the stationarity of the equilibrium average, the reversibility of the microscopic equations, and the fact that 
F(r a — r&) = F(||r a — 1*5 1|) is a central force. Secondly, in our work, the fluid is invariant under rotations by ip = |- 
for the hexagonal geometry and by (p = ^ for the square one. If we define the viscosity tensor as a linear operator f) 
acting on matrices A according to (^A)y = Tfo- /y A^\. Then our discrete symmetry can be written as: 



^(R -1 AiZ) = R _1 (j7A)R , 
for all matrices A, R being the rotation matrix 



(23) 
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cos if — sm if 
sin ip cos if 



(24) 



and f is equal to ^ or ^ respectively for the hexagonal or square systems. Thanks to this symmetry, the only 
nonvanishing elements are f)y : y = r\ji^j and rja^a — Vjj.jj- Furthermore, for i ^ j , k =/= I, 



— Vkl.kl j — 5 Viijj — T]kk,ll ■ (25) 

Hence, there are in fact only three independent elements: r\ xx . xx , rj xytXy , r) XXtyy . On the other hand, for an isotropic 
system, we can see that: 

V = Vxy^xy , (26) 

C = 2 (W* + ■ (27) 
The third element T) XX} y y is in fact a combination of the two other elements: 

(28) 



III. DESCRIPTION OF THE TWO-HARD-DISK MODEL 



In the present work, we apply our method to a simple model which we describe in the present section. The model 
is composed of two hard disks in elastic collisions on a torus. Bunimovich and Spohn have previously studied this 
model for a square geometry [10|. By periodicity, the system extends to a two-dimensional lattice made of infinitely 
many images of the two disks. For p. b. c. on a square domain, the infinite images form a square lattice, in which 
each cell contains two disks (see Fig. 




FIG. 1: The model of two hard disks: (a) in the hexagonal geometry and (b) in the square geometry. 



In the present work, we generalize this model to the hexagonal geometry (see Fig^a,). The possibility of such a 
model was pointed out in Ref. The images of each disk form now a triangular lattice. The two disks (the white 
and the black ones) have the same diameter a and mass m. They follow different trajectories. All the black disks 
move together and all the white ones move together. The system is periodic and the dynamics of the disks can be 
reduced to the dynamics in the unit cell or torus. 



A. Hexagonal geometry 



Let us first introduce some parameters of the system. L is the distance between the centers of two neighboring 
cells. It also corresponds to the distance between two opposite boundaries of a cell. 



FIG. 2: Basis vector (e and e'), position vector r a of particle a in the cell and the position vector v ala ; > in the lattice. 



By a linear combination of two vectors: 



fe = (L,0), 

we can spot all the cells of the lattice and then localize the center of a disk thanks to: 

r a l a i' a = r a + l a e + l' a e' , for a = 1,2, (30) 

where l a and l' a are integer, and r a is the position vector of the disk a with respect to the center of the cell. Therefore, 
the distance between the two disks is expressed by 



II ri h Vi - r 2 l2 ,i || = || r^-_r2 +(h - h) e + (l[ - 1' 2 ) e' || , (31) 

r 

where r = ri — r 2 is the relative position between both disks. By the minimum image convention, the relative distance 
||r|| should take the smallest value among the infinitely many possible values. Of course, this distance has to be 
greater than or equal to the disk diameter (||r|| = ||ri — r 2 | > er). As we have a hard-disk potential, the disks move 
in a free motion between each collision. Therefore, the equations of motion are written as: 





- El 


~dt 


m 


' dr 2 


P_2 


~~dt 


m 



6(t-t s ) , 



(32) 




(33) 

where pi and p 2 are the momenta of the two disks, Fi and F 2 being the forces applied respectively to the disks 1 
and 2. These forces equal zero when ||ri — r 2 || > a and are infinitely repulsive when ||ri — r 2 || = a. t s denotes the 
time of the jump to satisfy the minimum image convention. 
At this stage, we can do the following change of variables: 

C r = ri - r 2 , 

R= r L +r 2 (34) 
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Pi P2 

P - > (35) 

{ P = pi + p 2 . 

If we introduce the reduced mass \i = we can write: 

V% = P + I> Ar(s) =/iv + 5>Ar( s > S(t-t s ) , (36) 

S S 

|=F = Fl = -F 2 , (37) 

where v is the relative velocity and Ar( s ) = Ar^' — Ar^ . Here we suppose that we are in the reference frame of the 
mass center (that is P = 0). Accordingly, the energy of the system is reduced to: 



P 2 

£=£-. (38) 
2 /i 

The interest of this change of variables is to reduce the number of variables. Indeed, the only variables that remain 
are the relative position and velocity [r = [x, y) and v = (v x , v v )\. We can associate a fictitious pointlike particle with 
these variables, which moves in a reduced system, known as the periodic Sinai billiard (see Fig. |3J). 



d 




.twodisk: model. periodic Sinai billiard 



FIG. 3: The model of two hard disks in the hexagonal geometry is reduced to the periodic Sinai billiard thanks to a change of 
variables. 

The billiard is also a triangular lattice of hexagonal cells. The size d of these cells is equal to the size of the cells of 
the model itself (d = L). A hard disk is fixed on the center of each cell. Its radius is equal to the diameter a of the 
two moving disks. 

The basis vectors of this lattice are the same as those of the original dynamics (|32|l - l|33|l if we replace L by d, which 
gives us the possibility to spot a cell in the lattice thanks to the vector: 

r c = l c e + l' c e' , (39) 

where l c and l' c are integer. 

In the Sinai billiard, the system is described by a trajectory in a four-dimensional phase space which are the Cartesian 
coordinates (x,y,p x ,p y ), or the polar coordinates (x,y,pg,9). However, since the energy of the system is conserved, 
this space is reduced to the three-dimensional space of the variables (x, y, 6). Furthermore, in hard-ball systems, the 
topology of the trajectory is independent of the energy level. Therefore, we can study the system on an arbitrary 
energy level. This energy determines the temperature of the system and is equal to E — (d/2)(N — l)fcBT = k^T 
because we have only two degrees of freedom (d = 2, N = 2). Sinai and Bunimovich have demonstrated that the 
dynamics in such billiards is ergodic on each energy level [TH fl6| . 



9 



B. Square geometry 

The case of the square geometry is similar to the hexagonal one except that the basis vectors are here given by 

fe=(L,0), (4) 
\e' = (0,L), (4Uj 

where L is the length of a side of the square unit cell which contains two moving disks of diameter a. We perform 
the same change of variables to reduce the dynamics of two hard disks to the one of the fictitious pointlike particle of 
a Sinai billiard in a square unit cell. Here also, the size d of the cells of the Sinai billiard is the same as for the cells 
of the two hard disks model: d = L. 
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two- disk model periodic Sinai billiard 



FIG. 4: The model of two hard disks in the square geometry is reduced to the periodic Sinai billiard thanks to a change of 
variables. 



C. The different dynamical regimes of the model 

The physical quantity determining the size of the cell in our model is the density which corresponds to the number 
of disks per unit volume or, in our case, the number of disks per unit area. Each cell contains two disks. Therefore, 
the density is n — where V = ||e X e'|| is the area of a cell. In our study, we have chosen that the diameter of the 
moving disks is equal to the unity: a = 1. 

As a function of the density, we observe different dynamical regimes. At low density, the disks are able to move 
in the whole lattice so that the disks are not localized in bounded phase-space regions. In this case, the billiard may 
have a finite or an infinite horizon depending on the geometry and on the density. In the opposite, at high density, the 
disks are so close to each other that they cannot travel across the system and we refer to this regime as the localized 
regime. The critical density between the nonlocalized and localized regimes corresponds to the situation where both 
disks have a double contact with each other in the configuration shown in Fig. |SJ 




FIG. 5: Hexagonal system at the critical density n cr . 
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1. Hexagonal geometry 

In the hexagonal geometry, the area of the system is V — ||e x e'|| = ^-L 2 and the critical density is equal to: 

n cr = - 0.5774 , (41) 
even though the maximum density (the close-packing density) is: 

«max = 0.7698 . (42) 

y 

At the close-packing density, the system forms a triangular crystal. 

In the Sinai billiard, it is well known that there exists different kinds of regimes according to the dynamics of the 
particles. As a function of the density n, we observe three regimes: 

• The infinite- horizon regime: At the low densities < n < the particles can move in free flight over 
arbitrarily large distances. In this regime, the self-diffusion coefficient is infinite. (See Fig. 0) 




FIG. 6: Fypical configuration of the system in the infinite-horizon regime. 



• The finite- horizon regime: For the intermediate densities < n < n cr , the free flights between the collisions 
are always bounded by a finite distance of the order of the interdisk distance d. Therefore, the horizon is finite 
and the self-diffusion coefficient is positive and finite. (See Fig. 0) 




FIG. 7: Typical configuration of the system in the finite-horizon regime. 



• The localized regime: At the highest densities n cr < n < n max , the images of the disk overlap each other in the 
billiard so that the relative motion of the particles is localized in bounded regions. Therefore, the self-diffusion 
coefficient vanishes. (See Fig. |SJ) 

We notice that Figs. El [7J and [S] are not depicted at the same scale since the disk diameter is fixed to unity (<r = 1) 
and it is the interdisk distance d that varies. 

The infinite- and finite-horizon regimes extend over the densities < n < n CT . The localized regime corresponds 
to the densities n cr < n < n max . The following figure shows the different regimes in the hexagonal geometry. The 
remarkable feature of the hexagonal geometry is that there exists a finite-horizon regime which is not localized, in 
contrast to the square geometry (see below). 
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FIG. 8: Typical configuration of the system in the localized regime. 
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FIG. 9: The different dynamical regimes and thermodynamic phases of the model in the hexagonal geometry versus the density 
n. 

2. Square geometry 

In the square geometry, the volume is V = ||e x e'|| = L 2 and the critical density is: 

ricr = 0.5 , (43) 

which is the density of the transition between the infinite-horizon and the localized regimes. The close-packing density 
is equal to: 

«max = 1 • (44) 

In Fig. 1101 we have depicted the different regimes in the square geometry. In the square geometry, there also exist 
nonlocalized and localized regimes, but the horizon is always infinite in the nonlocalized regime. Therefore, it is only 
in the localized regime that the horizon is finite in the square geometry. This is an important difference with respect 
to the hexagonal geometry. 

IV. PROPERTIES OF THE MODEL 
A. Mean free path 

The mean free path (I) is the average distance between two successive collisions. It is known that, in two-dimensional 
billiards, the mean free path is related to the area A of the billiard and its perimeter C according to |l7j | 

(I) = ^ ■ (45) 



In the different regimes, the mean free path is given by 
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FIG. 10: The different dynamical regimes and thermodynamic phases of the model in the square geometry versus the density 
n. 
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FIG. 11: Theoretical (continuous line) and numerical (dots) values of the mean free path versus the density n in the hexagonal 
geometry. 



We show in Figs. 1111 and IT21 the excellent agreement between the above expressions and the values obtained by 
numerical simulations. The break observed in Figs. 1111 and IT21 between the nonlocalized and localized regimes can be 
explained thanks to Eq. Q45[l. Indeed, at the critical density n cr , the disks form a horn. Above criticality, the horn 
becomes a corner with a finite angle so that the perimeter C decreases very fast. But, on the other hand, the area A 
remains relatively constant. Therefore the ratio ^ increases with n until this effect disappears. At higher densities, 
the mean free path decreases again. 
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FIG. 12: Theoretical (continuous line) and numerical (dots) values of the mean free path versus the density n in the square 
geometry. 



B. Pressure and the different phases of the model 

The hydrostatic pressure allows us to interpret the different regimes in terms of thermodynamic phases. The 
pressure can be calculated in terms of the time average of the Helfand moment as shown in Appendix El In the 
two-disk model with N = 2 and d = 2, the pressure is given by 

PV = k B T + R , (46) 

where the rest can be calculated according to Eq. (|C8|) as 



R 4 (At-U ' (4?) 

where (At c _i c ) is the mean intercollisional time. If we denote by the angle between the velocity at collision and 
the normal to the disk of the Sinai billiard, the average in the numerator becomes 

{Ap[ c) • r$) = m v a (cos </> (c) ) , (48) 

a being the diameter of the disks. In the case the total momentum vanishes, the velocity v of the trajectory in the 
billiard is related to the relative momentum p, the energy, and the temperature by 

E = kBT= 2^ = m = — = — > (49) 
so that v = 2p/m. At collision, sin</()( c ' is uniformly distributed in the interval [— 1, +1] so that 

(cos^ c >> = £. (50) 

On the other hand, the mean intercollisional time of the billiard is related to the mean free path (/) and the speed 
v = 1 1 v 1 1 by 

(At c _ llC ) = ^ . (51) 
v 

Gathering the results, we obtain the rest as 

iramv 2 tt a , 

R = ^w = W) kBT - (52) 
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Accordingly, the hydrostatic pressure of the model is given by 



PV = k * T ( l + W))= kBT ( l + ^) ■ (53) 

In our work, we introduce the reduced pressure defined as 

P*= /3p}L = ™. = i + ^ = i + ^ (54) 

P N {N-\)k B T 4(/> AA [ ' 

In Figs. 1131 and ITU the reduced pressure is depicted as a function of the density and we observe the manifestation 
of a phase transition around the critical density. The hard-ball systems are known to present a fluid-solid phase 
transition that we here already observe in the two-disk model. 

At low density, the fictitious particle of the Sinai billiard can diffuse in the whole lattice. This means that the 
two disks move over arbitrarily large distances one with respect to the other, which is a feature of a fluid phase. In 
contrast, at high density, the fictitious particle is trapped between three (or four) disks and its motion is reminiscent 
of the vibration of atoms in a solid. Of course, it is not really a vibration since the disks bounce in a chaotic motion 
because of the elastic collisions whereas, in a solid, the atoms have quasi-harmonic oscillations around their equilibrium 
position. Nevertheless, we are in the presence of a solid phase because the translational invariance is broken. Indeed, 
the motion is no longer ergodic because the motion now is confined into one among several phase-space domains of 
the energy shell. 

A phase transition occurs between the fluid and solid phases. At the critical density n cr , the pressure has a 
maximum. Above ri cr , the pressure decreases, reaches a minimum at a value n' cr > n cr , before increasing again. For 
n CT < n < n' cx , the compressibility would be negative so that this state would be unstable from a thermodynamic 
viewpoint. This suggests a Maxwell construction to determine a fluid-solid coexistence in the interval of densities 
n-p < n < ng with rip < « cr and n' cr < ng. The values which would delimit this small coexistence interval in a 
thermodynamic interpretation of the transition would be given by 

• hexagonal geometry: np = 0.57 ±0.01, 

n s = 0.60 ±0.01 , (55) 

and 

• square geometry: np — 0.49 ± 0.01 , 

n s = 0.55 ±0.01 , (56) 

(see Figs. l9l and ITU)) . In the square geometry, the horizon is infinite in the fluid phase. In the hexagonal geometry, 
the horizon may also be finite in the fluid phase, which leads to finite viscosity coefficients in the fluid phase of this 
model as shown in the following. 




FIG. 13: Theoretical (continuous line) and numerical (dots) values of the reduced pressure P* versus the density n in the 
hexagonal geometry. 
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FIG. 14: Theoretical (continuous line) and numerical (dots) values of the reduced pressure P* versus the density n in the 
square geometry. 

V. VISCOSITY IN THE TWO-HARD-DISK MODEL 

A. The Helfand moment in the two-hard-disk model 

In our model defined with Eqs. 13411 and (|35|l and with a vanishing total momentum P = 0, the forces obey 
Fi = — F2 = F and the microscopic current can be written in relative coordinates as 



PlPj 



(57) 



where r is the smallest distance between the disks 1 and 2. Following the minimum image convention the position 
vector presents discontinuities because of the passages of the relative position through a boundary, after which it is 
reinjected into the cell at the opposite boundary. We denote the vectors normal to the boundaries of the unit cell by 



hexagonal geometry : < 



( C1 


= a , 




c 2 


= a 




c 3 


= b, 




c 4 


= b 


; 


c 5 


= b 


a . 


. c 6 


= a 


b. 



(58) 



and 



square geometry : 



= b 



(59) 



In order to satisfy the minimum image convention, the relative position undergoes jumps by vectors which are the 
vectors normal to the unit cell so that Ar( s ) = — c LJg where u s denotes the label of the boundary crossed by the 
particle at the s th passage at time t s . In these notations, Hamilton's equations take the form 



dr _ 2p 
dt m 
dp rp 
dt ~ r ■ 



Y,s c ^s s(t-t a ) , 



(60) 



In this periodic system, the expression for the Helfand moment is given by a reasoning similar to the one leading 
to Eq. (UHJ. We obtain: 



Gij (t) = Pi (t) r,(t) + Y,Mts) c ^s 3 0(t-t s ) 



(61) 
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Finally, the viscosity coefficients have the expressions: 

V13M = t lim o —= I / PiiQ <W? Pk(t s >) c Ub ,i )-(X>(*.) Cu,^ ) ( 51 C «V' / J • ( 62 ) 

~* 00 y \t«<t t a '« / \t s <t I \t B ,<t 1 1 

Let us remark that the terms Pi (t) rj (t) do not appear in this relation because they do not contribute to the viscosity 
coefficients. Indeed, the relative position r(t) and momentum p(t) remain bounded in the course of time and their 
contribution disappears in the limit t — *■ 00. 

In the following, the numerical results arc presented in terms of a reduced viscosity tensor which is defined by 

7 1ij,kl = ,-, / j \ rp ■ (63) 



B. Hexagonal geometry 

In the hexagonal geometry the fourth-order tensor of viscosity is isotropic. Indeed, since the system is invariant 
under rotations by |, we obtain the relation rj xx , yy = r} xx ,xx — 2 rjxy,xy which implies the full rotation invariance of 
the viscosity tensor. We depict in Figs. 1151 and the results obtained for the reduced viscosities (77*, £*) and the 
relation is checked in Fig. IT71 




FIG. 15: Shear viscosity coefficient rf versus the density in the hexagonal geometry. The part in dashed line corresponds to the 
density in which the coefficient would not exist in the limit t — > 00 because the horizon is infinite. The long dashed vertical lines 
separate the different regimes: on the left-hand side, the horizon- infinite regime (fluid phase); at the center, the horizon-finite 
regime (fluid phase); and on the right-hand side, localized regime (solid phase). 

In the infinite-horizon regime, the trajectory can present arbitrarily large displacements in the system without 
undergoing any collision. Accordingly, the variance of the Helfand moment G yx increases faster than linearly as 
t \ogt, which implies an infinite viscosity coefficient after averaging over an infinite time interval. However, the factor 
\ogt generates a so weak growth that it does not manifest itself much over the finite time of the simulation. This is 
the reason why we obtain finite values for the viscosity coefficients in Figs. H5l and Hrjl However, these values are only 
indicative since they should be infinite, strictly speaking. 

On the other hand, in the finite-horizon regime, the variance of the Helfand moment has a strictly linear increase 
in time and the viscosity coefficients are finite and positive. This is the result of a central-limit theorem which holds 
in the finite-horizon regime of the hexagonal geometry, as can be proved by considerations similar to those developed 
by Bunimovich and Spohn ^(| . We observe in Fig. that, the viscosity has a diverging singularity at the critical 

density (n cr = ^) which corresponds to the fluid-solid phase transition. We shall explain below the origin of this 
singularity. 
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FIG. 16: Bulk viscosity coefficient f* versus the density in the hexagonal geometry. 
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FIG. 17: Tensor element r]% y>xy of shear viscosity versus the density in the hexagonal geometry. The dots represent the results 
of the relation (1281 : Vxy.xy — \ {jfxx.xx ~ Vxx,yy)- The continuous line corresponds to the data of Fig. 1151 

Finally, in the localized regime corresponding to the solid phase, the viscosity is finite and positive, and decreases 
when the density increases until the maximum density. 



C. Square geometry 

In the square geometry, the fourth-order viscosity tensor is not isotropic. Indeed, the tensor is transformed by the 
matrix ijy (ip) of rotation by an angle <p into 



Vij.kli'P) = Rii'(<P) Rjj'iv) Rkk>{<f>) Rll'if) Vi']',k'l'(0) 



(64) 



For example, if (p = we have: 



?]xx,xx (f) = 2 \J]xx,xx (0) + Vxx,yy (0) + * Vxy,xy (0)] , 

Vxy,xy (f) = (0)] , 

?lxx,yy 

(f) = 

2 \J]xx,xx (0) + Vxx.yy (o)]- Vxy,xy (0). 



(65) 
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Since the system is not isotropic, one more viscosity coefficient is required beside the shear and bulk viscosities. 
Therefore, we have to evaluate the three independent tensor elements T) xx<xx , T) X y lXy , Tj xx>yy which are depicted in Figs. 
El and El with respect to two different axis frames: in the first one the axes are parallel to the sides of the square 
(ip = 0) and, in the second one, they form an angle of 45 degrees with respect to the lattice (ip = £). Figure HD1 shows 
that the relations (|o5l) are well satisfied between the elements of the viscosity tensor. 



(a) (b) (c) 




density density density 

FIG. 18: Square geometry: The three independent tensor elements (a) rj XXiXX , (b) r)% x ^ yy , (c) TjZ v , xv f° r P — 0. 




density density density 

FIG. 19: Square geometry: The three independent tensor elements (a) f] XXiXX , (b) T] xx>y y, (c) Vi%y, X y for ip =_§■ The continuous 
line corresponds to the results obtained numerically and the dots to the values obtained by the relations 1651 . 

In the square geometry, Bunimovich and Spohn have proved a central-limit theorem for viscosity in the localized 
regime which coincides with the solid phase above the critical density lOj. In this range of density, the viscosity 
coefficient is thus guaranteed to be positive and finite. 

In the fluid phase, the horizon is infinite and the viscosity is infinite because of a growth as ilogi of the variance 
of the Helfand moment for a reason similar as in the hexagonal geometry. In our numerical simulation over a finite 
time interval, the viscosity takes finite values because the logarithmic growth is very slow. 

An important difference with respect to the hexagonal geometry is the absence of a singularity of the viscosity 
coefficient T)* yxy (Q) at the phase transition in the square geometry. However, such a singularity still appears in the 
square geometry in the coefficients r] xxxx (0), r] xx yy (0), and vl v , xy (j )■ 

Moreover, in the solid phase, the coefficient rf xy xy (0) increases with the density, as explained here below. 

D. Explanation of the numerical observations 

1. Solid phase 

The behavior of the viscosity tensor is clearly different in the two geometries. In this section, we explain these 
differences by comparing the topology of the trajectories in both geometries, since these trajectories form the basis of 
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the evolution of the Helfand moment. More precisely, we will compare the behavior of rf xy between the hexagonal 
and square geometries for ip = 0. This viscosity coefficient is given by 

Vxy/xy — Vyx,yx ~ " i ~ ' (Gyx(t)) = , (66) 

Gyx ^ ^^s^y(^s) ^lj s x j 

where v y (t s ) is the y-component of the velocity at the time t a of the jump. 

When the density tends to the closed-packing density, the accessible domain of the particles tends to a perfect 
triangle in the hexagonal geometry. On the other hand, in the square geometry, it tends to a perfect square. This 
difference is at the origin of the different behaviors of the r)% yxy m both lattices. 



y 




X 



FIG. 20: Part of a typical trajectory in the square geometry when the density tends to the closed-packing density. 

First, let us consider the case of the square geometry. In Fig. [201 we depict a typical trajectory of the fictitious 
particle moving in the Sinai billiard. We observe that this trajectory presents a regular motion between two opposite 
"walls" (these walls are made of parts of the fixed hard disks in the billiard) . At the limit where the billiard is a perfect 
square, the trajectories will bounce back and forth in a regular motion. Indeed the square billiard is an integrable 
system. 




FIG. 21: Geometry and notation for the boundaries in the case of the square geometry at high density. 

As we have seen before, the evolution of the Helfand moment along the trajectories is determined by the passages 
through the boundaries (see Fig. I21[l . Both horizontal boundaries (3 and 4) do not contribute to the evolution of G yx 
since the x-component of the normal vectors to these boundaries equals zero. Therefore, only the passages through 
the vertical boundaries contribute to the Helfand moment in the square geometry. 

To understand the behavior of the Helfand moment, let us take a small part of the typical trajectory drawn in Fig. 
1201 (see Fig. I22JI . First, let us consider the part denoted by the letter a in Fig. [221 This one crosses the boundary in 
the direction 1 — > 2, which means that c LOsX is positive (since c% x = 4). On the other hand, the y-component of the 
velocity, v y , is also positive. Therefore, the contribution of the small part a to the evolution of G yx is positive. 

Now, let us take the part of the trajectory denoted b in Fig. [221 In this case, the particle crosses the boundary 
in the direction 2 — * 1 and C2x is negative. Since v y is also negative, the product of these two quantities is positive, 
and so at each successive crossings of the boundary 1 — 2. Consequently, we obtain a sum of positive terms and the 
Helfand moment quickly increases along a trajectory as the one of Fig. |2~U1 

However, the square is not perfect and the walls are still slightly convex. Therefore, after a certain time, the 
trajectory shown in Fig. I2UI goes into a transient regime shown on the left-hand side of Fig. 1231 before another regime 
in which the particle collides most often the two other walls (see the right-hand side of Fig. I23fl . 
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FIG. 22: Part of a typical trajectory in the square geometry at high density. 



M 00 




FIG. 23: Square geometry at high density: The trajectory is depicted (a) during a transient regime before (b) another regime 
with most bounces on the two other opposite walls. 



With the same reasoning as before, we conclude that the contributions are negative in this new regime and the 
Helfand moment decreases during a long time interval. 

The evolution of the Helfand moment along the whole trajectory is depicted in Fig. |2] where we observe the 
succession of the three types of regimes which we have described here above. We notice that the nearly constant part 
corresponds to the transient regime. 

The larger is the density the more perfect is the square and the longer the trajectory remains in a particular regime. 
Therefore, the Helfand moment can have larger and larger variations, which implies an increase of the coefficient 
T]xy iXy (0) of shear viscosity with density. 

In the hexagonal geometry (see Fig. I25|) . the trajectories present another behavior. We show in Fig. [23] a typical 
trajectory in this geometry with a density larger than the critical density. We observe that the trajectory visits the 
whole billiard in different directions and therefore goes into very different velocities. Accordingly, the particle crosses 
the boundaries with random values of its velocity in contrast to its behavior in the square geometry. Consequently, 
the quantity can be positive at a particular crossing and negative at the next one. Hence the Helfand moment 
cannot increase or decrease over long periods as in the square geometry (see Fig. This explains qualitatively 

why, in the solid phase, the coefficient rj* y xy (0) = rf is much smaller in the hexagonal geometry than in the square 
one. 

In the square geometry with ip = j, the same arguments as in the hexagonal case explain the decrease of r/* y xy (j) 
at high density. By the relations between the different elements of the viscosity tensor, we can also understand the 
behavior of the other elements in both geometries. 



2. Fluid-solid phase transition 

In both the hexagonal and square geometries, the two-disk model presents a phase transition. This transition is 
reminiscent of the fluid-solid phase transition in the many-disk system where the viscosity coefficient is also singular. 
In this regard, the two-disk model can contribute to the understanding of the changes in the transport properties 



21 



50QQ 



5000 



-10000 



HXKX) 2<XXX> MM® 4<KXX> 

collision number 



FIG. 24: Evolution of the Helfand moment along a typical trajectory in the square geometry at high density. 
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FIG. 25: Geometry and notation for the boundaries in the case of the hexagonal geometry at high density. 

across the fluid-solid phase transition. 

We first explain why rf xy xy presents a diverging singularity at the critical density in the hexagonal geometry and 
not in the square geometry for ip = 0. Here again, we compare the topology of the trajectories in both geometries 
and the way in which the Helfand moment evolves along these trajectories. At densities close to the critical density, 
both geometries present what we call traps. 

Figure 1271 shows an example of a trap. These traps are particular regions of the billiard where the particle can 
remain during a long time interval. Figure depicts typical examples of a particle moving in such traps. When the 
particle travels out of the traps, the Helfand moment does not increase quickly in both geometries. Therefore, it is 
the presence of the traps which is at the origin of the difference between both geometries. 

In the square geometry, the traps do not influence the evolution of G yx . Indeed, as we have already mentioned here 
above, the passages through the horizontal boundary 3 — 4 do not contribute since cs x = c^ x = (see Fig. |^for the 
definitions of the boundaries in the square geometry). Therefore, the horizontal traps around these boundaries do 
not contribute. There remains the vertical traps. When a particle bounces for a long time in one of these traps, c\ x 
and ci x are not vanishing, but the velocity v y is almost equal to zero so that the vertical traps does not contribute 
much either. This implies that both kinds of traps contribute very slightly to the evolution of the Helfand moment. 
To conclude the Helfand moment diffuse in the same way as for the other densities and the coefficient rj* xy (0) does 
not present any divergence at the fluid-solid transition in the square geometry. 



On the other hand, in the hexagonal geometry, the traps along the boundaries making an angle of 30° with respect to 
the horizontal are very important for the evolution of G yx , whereas the vertical traps do not participate significantly. 
Figure [20 shows a typical diffusion of the Helfand moment. We observe in Fig. [20 the presence of jumps which 
correspond to the passages in the traps like the one drawn on the left-hand side of Fig. [20 Because of these jumps, 
the Helfand moment quickly diffuses. Furthermore, the importance of these traps in the hexagonal geometry can also 
be understood by comparing the behavior of the Helfand moment as a function of time at densities below and above 
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FIG. 26: Hexagonal geometry at high density: (a) Part of a typical trajectory when the density tends to the closed-packing 
density, (b) Evolution of the Helfand moment along this typical trajectory. 




FIG. 27: Example of traps in which the particles can enter and remain a long time. 



the critical one n cr . 

We illustrate this point in Fig. 1301 where we observe that there are no more jumps above the critical density. 
Therefore, G yx does not vary much contrary to the case of densities just below n cr . Above criticality, the size of the 
traps decreases so quickly that the contribution of these traps decreases and, thus, the viscosity coefficient rj* y xy = if 
also decreases. By these arguments, we have an explanation for the diverging singularity of the shear viscosity at the 
phase transition in the hexagonal geometry. 

This results show that, at a fluid-solid phase transition the viscosity coefficients may depend sensitively on the 
geometry of the lattice of the solid phase in formation. 



E. Viscosity by the method of Alder et al. 

We have also verified numerically that our method of calculation of the viscosity based on the Helfand moment 
(|16fl gives the same values as the method of Alder et al. based on the expression (|21|l In the two-disk system, this 
expression reduces to: 



(67) 



As shown in Fig. |^for the shear viscosity in the hexagonal geometry, there is an excellent agreement between the 
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FIG. 28: Particle trapped between two disks very close to each other in the hexagonal geometry. The line joining their centers 
either (a) forms an angle with the horizontal or (b) is horizontal. 
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FIG. 29: Helfand moment in the hexagonal geometry evaluated along a particular trajectory at a density tending to the critical 
density. 

values obtained by both methods, which confirms the exact equivalence of both methods. 



VI. VISCOSITY IN SYSTEMS OF N HARD DISKS 

In the present section, we apply our Helfand-moment method to systems of N hard disks. Our purpose is to show 
that the values of the shear viscosity obtained for the two-hard-disk model are in good agreement with the values for 
larger systems, as well as with Enskog's theory. 

Figure |23 depicts the shear viscosity of systems containing from N — 2 up to N = 40 hard disks. For TV = 2, 
we consider here the hexagonal geometry. For the systems with N = 4-40 disks, the time evolution is simulated by 
molecular dynamics with periodic boundary conditions in the square geometry. The viscosity is calculated by the 
Helfand-moment method based on Eq. (|21(1 . 

We observe that, at low densities, the numerical values are in very good agreement between themselves. At higher 
densities, differences appear because the fluid-solid transition shifts toward higher densities as the number of disks 
increases. For N = 2 disks in the hexagonal geometry, the fluid-solid transition occurs in the interval n — 0.57-0.60, 
while it occurs in the interval n — 0.87-0.90 for N = 40. The sharp singularity of viscosity for N = 2 in the hexagonal 
geometry is specific to the geometrical constraints of a two-degree-of-freedom system, as explained in the previous 
section. Nevertheless, we notice that the decrease of the shear viscosity just above the fluid-solid transition is also the 
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FIG. 30: Comparison of the evolution of the Helfand moment for two different densities separated by the critical density in the 
hexagonal geometry. 




FIG. 31: Shear viscosity rf in the hexagonal geometry calculated by our Helfand moment l|16|l (continuous line) and the one 
of Alder el al. (dots). 



feature of the large system with TV = 40 disks. 

Furthermore, the results of our Helfand-moment method are compared with Enskog's theory. For a fluid of hard 
disks of mass m and diameter a, Enskog's theory predicts that the shear viscosity is given by [3 



V = Vo ( y + 2y + 3A9l6Y y' ) , (68) 



where 



1.022 mk B T 

Vo = ^r—\ — - — ■ (69) 
2 a V 7r 

is the Boltzmann value of the shear viscosity, Y is the Enskog factor entering the equation of state as follows: 



P = nk B T(l + 2yY) 



(70) 
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FIG. 32: Shear viscosity r\ versus particle density n in fluids at temperature T = 1 with N — 2,4, 8, 12, 16, 20, 40 hard disks of 
unit mass and diameter. The solid line is Enskog's value l|68fl . For N — 2, the data are the same as in Fig. 1151 except that we 
here plot r\ — 2rj* instead of rf as in Fig. 1151 



and y = ira 2 n/4. For the hard-disk fluid, a good approximation of the Enskog factor is given below the fluid-solid 
transition by |19j : 

1 Z- v 

Y = (71) 

It is known that the Enskog approximation is not good around the fluid-solid transition and at very high densities. 

A remark is here in order. It is known [20| that the viscosity coefficient of the infinite hard-disk fluid is diverging 
because of long-time tails. However, this divergence is only logarithmic and does not manifest itself in numerical 
calculations before extremely long times. This explains why the long-time tails do not spoil the agreement between 
the numerical values and Enskog's theory. 

We see in Fig. [22 the good agreement between Enskog's theory and the numerical values of our Hclfand-momcnt 
method at low densities showing the consistency of our results. 



VII. CONCLUSIONS 



In this paper, we propose a new expression for the Helfand moment associated with viscosity in molecular dynamics 
with periodic boundary conditions. This new Helfand moment takes into account the minimum image convention 
at the basis of molecular-dynamics simulations with periodic boundary conditions. In order to satisfy the minimum 
image convention, the position coordinates of the particles undergo jumps. These jumps modify both the equations 
of motion and the Helfand moment which is given by the time integral of the microscopic current entering the Green- 
Kubo formula. As a consequence, the viscosity tensor calculated with our Helfand moment is equivalent to the one 
based on the Green-Kubo formula, as proved in Appendix [5] I n the case of hard-ball systems, we also prove in 
Subsec. Ill Dl that our method is equivalent to the method by Alder et al. Q. Moreover, we show in Appendix [Q that 
the hydrostatic pressure can also be calculated thanks to the Helfand moment we propose for molecular dynamics 
with periodic boundary conditions. Our new Helfand moment and our proofs bring a solution to the ambiguities and 
problems reported by Erpenbeck about the definition of a Helfand moment in a molecular dynamics with periodic 
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boundary conditions. We think that the Helfand-moment method can be very useful for the numerical calculation of 
viscosity because this method has the advantage of being numerically robust. 

We have applied our Helfand-moment method to the numerical calculation of viscosity in systems of hard disks. 

In Sec. viscosity has been studied in detail in a simple model composed of two hard disks in elastic collision. 
This model has already been investigated in the square geometry by Bunimovich and Spohn |Ioj|. In the present 
paper, we generalize this model to the hexagonal geometry. First, we show that the fourth-order viscosity tensor is 
isotropic in the hexagonal geometry although it is not in the square geometry. Secondly, we show the viscosity can 
be positive and finite in the fluid phase of the hexagonal geometry, although it is always infinite in the fluid phase of 
the square geometry. The reason is that the horizon of the Sinai billiard driving the dynamics of the two-disk model 
is always infinite in the fluid phase of the square geometry although there is a regime with a finite horizon in the fluid 
phase of the hexagonal geometry. In an infinite-horizon regime, the viscosity becomes infinite so that, from a physical 
point of view, the proof of the existence of a positive and finite viscosity coefficient strictly holds in the hexagonal 
two-disk model. In the solid phase, the transport coefficients acquire a different meaning because the spontaneous 
breaking of translational invariance modifies the structure of the hydrodynamic modes and the viscosity coefficient 
should be reinterpreted in terms of the damping coefficients of the transverse sound modes and of the diffusive modes 
pll I22I |23| | . We hope to report on this question in a future publication. 

The two-disk model presents a phase transition between a fluid and a solid phase. This transition is reminiscent 
of the fluid-solid transition in the system composed of many disks. Indeed, the transition manifests itself in the 
hydrostatic pressure in a very similar way as in the many-particle system. The hydrostatic pressure can be directly 
related to the mean free path in the two-disk model and we can thus explain the manifestation of the transition on the 
pressure in terms of the behavior of the mean free path near the transition. In this simple model, the transition can 
be understood as a geometrical property of the dynamical system. Indeed, the trajectories are unbounded in the fluid 
phase albeit there remain localized in bounded domains in the solid phase where ergodicity is broken. The fluid-solid 
transition also manifests itself as a diverging singularity in the viscosity in the two-disk model. We have here shown 
that this singularity in the viscosity versus the density may depend sensitively on the geometry of the lattice of the 
solid phase in formation. 

In Sec. IVII we have extended the calculation of shear viscosity to systems with many disks. The remarkable result is 
that the two-disk systems already gives the shear viscosity in quantitative agreement with its values in larger systems, 
as well as with Enskog's theory at moderate densities. 

In a companion paper, we report a study of viscosity by the escape-rate method [24| . In this other work, we use 
the Helfand moment which we have introduced in the present paper. 
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APPENDIX A: MICROSCOPIC DERIVATION OF THE VISCOSITY TENSOR 



In this Appendix, we provide a short microscopic derivation of the viscosity tensor. 

First, we need the balance equation for the local conservation of momentum. If we define the density of momentum 

as 

N 

9i(r) = ^2 Pat 5 (r-r a ) , (Al) 

0=1 

the balance equation is 

d t 9i + d 3 T i:i = , (A2) 
with dj = d/drj. The microscopic momentum current density is given by 

N N 1 

r . . = £ E^l S (r -r a ) + ±J2 *K r « - ^ / dX 6 t r - r ^ ' ( A 3) 

where r a h(A) is the parametric equation of a curve joining the particles a and b: r a (,(0) = and r a fc(l) = r a . 
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The microscopic current associated with viscosity is defined by integrating the momentum current density over the 
volume V: 



Jij = / Tij (r) dr 
Jv 



which is given by Eq. ©. We notice that the hydrostatic pressure is given at equilibrium by 



(Jij^eq — P V Sij , 



(A4) 



(A5) 



if second-order tensors are isotropic in the system of interest. 

We suppose that, at the initial time, the fluid is close to the equilibrium and described by the following nonequi- 
librium distribution: 



V(T) = V eq (T) 



1 + J g(r) • v(r) dr 



= Pcq(r) 



N 



l+/?^p a -v(r a ) 



where V cq is the equilibrium distribution and is a normalization constant such that 



(Pai Pbj) cq ^ &ij &ab 



In the microcanonical state, we have that 



0- 



1 N 



kv,T N - 1 



(A6) 



(A7) 



(A8) 



The aforementioned distribution describes a fluid with a macroscopic velocity field v(r) since the nonequilibrium 
average of the momentum density can easily be shown to be given by 



<g(r)). 



where 



Pcq v(r) 



N 



p cq = m — , 



(A9) 



(A10) 



is the mass density at equilibrium. 

The time evolution of the probability density (|A6|) is ruled by Liouville's operator given by the Poisson bracket 
with the Hamiltonian L = {H, •} or the pseudo-Liouville operator in the case of hard-ball dynamics. This operator 
has the effect of replacing the phase-space coordinates T by r(— t) 



V t - e Lt Vo = V cq (T) 



1 + e g{r) ■ v(r) dr 



^c q (r) 



N 



l+^p a H).v[r a H)] 



(All) 



Alternatively, we known that the time evolution of the momentum density is given by Eq. (|A2|) . In this case, the 
momentum density should be considered as an observable so that the solution of Eq. (|A2|I is 



so that 



is solution of the equation 



Integrating both sides over time we get 



g(r,t) =e- Zt g(r,0) , 
e £ *g(r) =g(r,-i) , 
d t 9i = dj nj . 

ft(r,-t) = fli(r,0)+ / dt'd jnj (t') 



(A12) 
(A13) 
(A14) 

(A15) 
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Close to equilibrium, we may consider the time evolution of deviations with respect to the equilibrium. We neglect 
terms which are quadratic in the deviations such as the velocity field itself. The time evolution of these deviations is 
obtained by considering the nonequilibrium average of the balance equation (|A2J) for the deviations: 

d t (%)noncq + 9j (STij) no neq = , (A16) 

with 

STij = Tij - (Ty) oq . (A17) 

The nonequilibrium average of the deviation of the momentum current density is given by 

(<Mr)) noncq = J 5 nj (r) V(T, t) dT = f3 J dr'(5 nj (r) g k (r' , -t)) eq v k (r') . (A18) 
We use Eq. (| Al 5|> to transform the average as 

(Sr l3 (r) g k (r',-t)) eq = (St zj (v) . 9fc (r',0)) oq + / dt' (6 Tij (r,0) d[ 5t m (r', t')) oq , (A19) 



where we have used the property that d[(T k i) cq = because the equilibrium state is spatially uniform. We notice 
that the first term in the right-hand side of Eq. (|A19|) vanishes because the equilibrium average of an odd power of 
particle momenta vanishes. After an integration by part over the velocity field, Eq. I|A18(1 becomes 

(<Mr)) nonoq = -13 J dr' J dt' (*ry(r,0) 8 m {v\ i')) cq d[v k {v') = - VljM div k (r) , (A20) 
where the identification with the viscosity tensor is carried out in the limit t — > oo by 

mjM S(r -r')=P dt' (<Mr,0) ^(r', <')>c q . (A21) 



Taking the double volume integral J v dr J v dr' of both sides of Eq. I|A21|) and dividing by the volume V, we obtain 
the viscosity tensor as 

f°° 

%',« = 77/ dt ( SJ ij(.°) SJki(t)) cq , (A22) 



VJ 

with 



SJij(t) = I dr ST tJ (r, t) = J l3 (t) - (J y ) eq , (A23) 

Q.E.D. 



APPENDIX B: PROOF OF THE EQUIVALENCE BETWEEN GREEN-KUBO AND 

EINSTEIN-HELFAND FORMULAS 

Our aim is here to deduce the Green-Kubo formula JJJ from the Einstein-Helfand formula (|14|l , proving the equiva- 
lence between both formulas under the condition that the Helfand moment is defined by Eq. I|15|l as the time integral 
of the microscopic current (JSJ) and the further condition that the time auto-correlation functions decrease fast enough. 

We start from the Einstein-Helfand formula (|14|) with 

SG Z] (t) = f SJij(T) dr , (Bl) 
Jo 

SJij being defined by Eq. I|A23|I and supposing for simplicity that SGij(0) = 0. Accordingly, we have successively 
from Eq. (O that 

miM - T lim JL (5G l3 (T)SG kl (T)) 
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lim Jj- I dt x [ dt 2 (SJ^ihjSMh)) 

T^oo Zl V J Q 



T r T-\t\/2 



lim ^ I dt I dr (5J l3 (0)5J kl (t)) 

T^oo Zl V J_ T J\t\/2 



-\-oo 



^ I <H ;<U, ( (UM. /,,(/!) 



-oo 

p_ r + °° 

V Jo 



dt (6^(0)6 J u (t)) , (B2) 



where we have performed the change of integration variables 



t = t 2 -h 

rr- — h+iS 

' 2 ' 



(B3) 



and supposed that 



1 



+T 



lim - / dt \t\ (6^(0)8 J kl (t)) = , (B4) 



T 



which requires that the time autocorrelation functions decrease faster than \t\ 1 e with e > 0. Q.E.D. 

APPENDIX C: PRESSURE AND HELFAND MOMENT 

The hydrostatic pressure at equilibrium is given as the mean value of the momentum current density, i.e., as the 
mean value of the same microscopic current entering the Green-Kubo relation: 



PijV = / (7y) eq dr = (Jy)eq . (CI) 
Jv 

The average over the thermodynamic equilibrium state can be replaced by a time average: 

1 I* 

PijV = (J tJ ) cq = lim - / dr Jij . (C2) 



(i 



t^oo t 

We can here introduce the Helfand moment to obtain the hydrostatic pressure from the Helfand moment as 



In the microcanonical equilibrium state we have that 

N - 1 

[Pai Paj)eq = m k B T — 5y . (C4) 

If we assume that the system is isotropic, Pij = P Sij and we obtain 

PV = {N-l)k B T + R, (C5) 

where the rest R provides the corrections to the law of perfect gases in dense systems. By using Eqs. (|16[) and 121JI . 
the virial can be computed alternatively by 

JV 

/ i 

R = 



^ £ F(r afc )-r Qb \ (C6) 

= ^^EEp^-^^-^) (C7) 

s a— 1 

= ^^£ A P£ c) -r$*(*-0 (C8) 
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where d is the dimension, r a t — r a — r^, t s are the times of jumps to satisfy the minimum image convention, while 
the last expression only holds for hard-ball systems, t c are the collision times, and r^j} = r a (t c ) — rb(t c ). 
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